Galton Board

WLJS Version adapted from [Galton Board

2025, Gabriel Lemieux](https://www.wolframcloud.com/obj/gabriel.lemieux/Published/GaltonBoard.nb)

(* Initialization *) ball = <| r -> 0.1, (* Radius of the ball *) n -> 100 (* Number of simultaneous balls *) |>; peg = <| r -> 0.25 (* Radius of the pegs *) |>; sim = <| \[CapitalDelta]t -> 0.05, (* Higher Value need to correct collision positions *) \[Epsilon] -> 0.4, (* 1.0 means all energy back in the ball *) g -> {0, -9.81}, (* You can tilt the board *) n -> 1000 (* Total Number of Balls *) |>; (* Initialize Position and Momentum Matrices *) randomBall := {RandomReal[{-0.05, 0.05}], RandomReal[{0, 1}] + 10}; ball[pos] = Table[randomBall, ball[n]]; ball[mom] = Table[0, ball[n], 2]; ball[bins] = {}; (* Peg position *) peg[pos] = Flatten[Table[{ii + 0.5 jj, jj}, {ii, -15, 10}, {jj, -1, 8}], 1]; (* Display function *) positions = ball[pos]; Graphics[{ {LightGray, EdgeForm[Gray], Disk[{#1, #2}, peg[r]]} & @@@ peg[pos], {PointSize[0.02], Opacity[0.8], Blue, Point[positions // Offload]}, AnimationFrameListener[positions // Offload, "Event"->"frameTick"] }, AspectRatio -> 0.5, PlotRange -> {10 {-1, 1}, 10 {0, 1}}, PlotRangeClipping -> True, ImageSize -> Large ] EventHandler["frameTick", Function[Null, If[Length[ball[pos]] > 0, (* Step *) (* Scheme: Backward Euler *) (* Next Step Momentum *) ball[mom] += Threaded[sim[\[CapitalDelta]t] sim[g]]; (* Next Position Prediction *) next = ball[pos] + sim[\[CapitalDelta]t] ball[mom]; (* Distance: Do not parallelize, it's already in the algo. *) dist = DistanceMatrix[next, peg[pos]]; (* Collisions between current and future step *) col = Parallelize @ MapIndexed[ {v, i} |-> { SelectFirst[v, (# < ball[r] + peg[r]) & -> "Index"], i[[1]] } /. ({Missing[_], _} -> Nothing), dist ]; (* Correct Momentum to include collisions *) If[Length[col] > 1, { up, val } = Transpose[ Parallelize @ MapApply[ {p, b} |-> { b, sim[\[Epsilon]] * Norm[ball[mom][[b]]] Normalize[next[[b]] - peg[pos][[p]]] }, col ] ]; momt = ball[mom]; momt[[up]] = val; ball[mom] = momt; ]; (* Update pos with corrected momentum *) ball[pos] += sim[\[CapitalDelta]t] ball[mom]; (* Filter through the balls, remove/reset and register those exiting *) exited = Select[ball[pos], #[[2]] < 0 & -> "Index"]; ingame = Select[ball[pos], #[[2]] > 0 & -> "Index"]; If[Length[exited] > 0, ball[bins] = Join[ball[bins], ball[pos][[exited, 1]]] ]; (* Remove them from the pos and mom matrices *) ball[pos] = ball[pos][[ingame]]; ball[mom] = ball[mom][[ingame]]; (* Insert new ball *) left = sim[n] - Length[ball[bins]] - Length[ball[pos]]; toAdd = Min[left, ball[n] - Length[ball[pos]]]; Do[ AppendTo[ball[pos], randomBall]; AppendTo[ball[mom], {0, 0}], toAdd ]; positions = ball[pos]; ]]];